home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / poly_fit.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  135 lines
  1. ; $Id: poly_fit.pro,v 1.6 1997/01/30 19:21:24 ali Exp $
  2.  
  3. FUNCTION POLY_FIT, X, Y, NDEGREE, YFIT, YBAND, SIGMA, CORRM, DOUBLE=double
  4. ;+
  5. ; NAME:
  6. ;    POLY_FIT
  7. ;
  8. ; PURPOSE:
  9. ;    Perform a least-square polynomial fit with optional error estimates.
  10. ;
  11. ;    This routine uses matrix inversion.  A newer version of this routine,
  12. ;    SVDFIT, uses Singular Value Decomposition.  The SVD technique is more
  13. ;    flexible, but slower.
  14. ;
  15. ;    Another version of this routine, POLYFITW, performs a weighted
  16. ;    least square fit.
  17. ;
  18. ; CATEGORY:
  19. ;    Curve fitting.
  20. ;
  21. ; CALLING SEQUENCE:
  22. ;    Result = POLY_FIT(X, Y, NDegree [,Yfit, Yband, Sigma, CORRM] )
  23. ;
  24. ; INPUTS:
  25. ;    X:    The independent variable vector.
  26. ;
  27. ;    Y:    The dependent variable vector, should be same length as x.
  28. ;
  29. ;     NDegree:    The degree of the polynomial to fit.
  30. ;
  31. ; OUTPUTS:
  32. ;    POLY_FIT returns a vector of coefficients with a length of NDegree+1.
  33. ;
  34. ; OPTIONAL OUTPUT PARAMETERS:
  35. ;    Yfit:    The vector of calculated Y's.  These values have an error 
  36. ;        of + or - Yband.
  37. ;
  38. ;    Yband:    Error estimate for each point = 1 sigma
  39. ;
  40. ;    Sigma:    The standard deviations of the returned coefficients.
  41. ;
  42. ;    Corrm:    Correlation matrix of the coefficients.
  43. ;
  44. ; Keyword Parameters:
  45. ;    DOUBLE = if set, force computations to be in double precision.
  46. ; COMMON BLOCKS:
  47. ;    None.
  48. ;
  49. ; SIDE EFFECTS:
  50. ;    None.
  51. ;
  52. ; MODIFICATION HISTORY:
  53. ;    Written by: George Lawrence, LASP, University of Colorado,
  54. ;        December, 1981.
  55. ;
  56. ;    Adapted to VAX IDL by: David Stern, Jan, 1982.
  57. ;       Modified:    GGS, RSI, March 1996
  58. ;                    Corrected a condition which explicitly converted all
  59. ;                    internal variables to single-precision float.
  60. ;                    Added support for double-precision inputs.
  61. ;                    Added a check for singular array inversion.
  62. ;             SVP, RSI, June 1996
  63. ;                     Changed A to Corrm to match IDL5.0 docs.
  64. ;
  65. ;-
  66.     ON_ERROR,2        ;RETURN TO CALLER IF ERROR
  67.  
  68.     N = N_ELEMENTS(X)     ;SIZE
  69.     IF N NE N_ELEMENTS(Y) THEN $
  70.       message,'X and Y must have same # of elements'
  71. ;
  72.     M = NDEGREE + 1L ;# OF ELEMENTS IN COEFF VEC.
  73. ;
  74.  
  75.   sx = size(x)
  76.   sy = size(y)
  77.   if n_elements(double) eq 0 then $    ;True if working in double prec
  78.     double = (sx[sx[0]+1] eq 5) or (sy[sy[0]+1] eq 5)
  79.  
  80.   if double then begin
  81.     CORRM = DBLARR(M,M)        ;Operands are double, COEFF MATRIX
  82.     B = DBLARR(M)        ;WILL CONTAIN SUM Y * X^J
  83.     Z = Replicate(1.0d0, N)
  84.         xx = x
  85.   endif else begin
  86.         CORRM = fltarr(M,M)         ;COEFF MATRIX
  87.         B = fltarr(M)           ;WILL CONTAIN SUM Y * X^J
  88.         Z = Replicate(1.0, N)
  89.         xx = float(x)
  90.   endelse
  91. ;
  92.     B[0] = TOTAL(Y, DOUBLE = double)
  93.     CORRM[0,0] = N
  94. ;
  95.     FOR P = 1,2*NDEGREE DO BEGIN ;POWER LOOP.
  96.         Z=Z*XX            ;Z IS NOW X^P
  97.                 ;B IS SUM Y*X\X^J
  98.         IF P LT M THEN B[P] = TOTAL(Y*Z)
  99.         SUM = TOTAL(Z)
  100.         FOR J= 0 > (P-NDEGREE), NDEGREE < P DO CORRM[J,P-J] = SUM
  101.       END            ;END OF P LOOP.
  102. ;
  103.     CORRM = INVERT(CORRM, status)    ;INVERT MATRIX.
  104.         if status ne 0 then message, "Singular matrix detected."
  105. ;
  106. ;            IF CORRM IS MULTIPLIED BY SIGMA SQUARED, IT IS THE
  107. ;            CORRELATION MATRIX.
  108. ;
  109.     C = b # CORRM    ;Get coefficients
  110. ;
  111.     IF (N_PARAMS(0) LE 3) THEN RETURN,C    ;EXIT IF NO ERROR ESTIMATES.
  112. ;
  113.     YFIT = REPLICATE(C[0], N)  ;Initial Yfit
  114.     FOR K = 1,NDEGREE DO YFIT = YFIT + C[K]*(XX^K) ;FORM YFIT.
  115. ;
  116.     IF (N_PARAMS(0) LE 4) THEN RETURN,C    ;EXIT IF NO ERROR ESTIMATES.
  117. ;
  118.     IF N GT M THEN $ ;COMPUTE SIGMA
  119.         SIGMA = TOTAL((YFIT-Y)^2) / (N-M) ELSE    SIGMA = 0.
  120. ;
  121.     CORRM=CORRM* SIGMA        ;GET CORRELATION MATRIX
  122. ;
  123.     SIGMA = SQRT(SIGMA)
  124.     YBAND = FLTARR(N)+ CORRM[0,0]    ;SQUARED ERROR ESTIMATES
  125. ;
  126.     FOR P = 1,2*NDEGREE DO BEGIN
  127.       Z = XX ^ P
  128.       SUM = 0.
  129.       FOR J=0 > (P - NDEGREE), NDEGREE < P DO SUM = SUM + CORRM[J,P-J]
  130.       YBAND = YBAND + SUM * Z ;ADD IN ERRORS.
  131.     END        ;END OF P LOOP
  132.     YBAND = SQRT(ABS(YBAND))    ;ERROR ESTIMATES
  133.     RETURN,C
  134. END
  135.